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Abstract 



Disease control is of paramount importance in public health with infectious disease extinction as the 
ultimate goal. Although diseases may go extinct due to random loss of effective contacts where the 
infection is transmitted to new susceptible individuals, the time to extinction in the absence of control may 
be prohibitively long. Thus intervention controls, such as vaccination of susceptible individuals and/or 
treatment of infectives, are typically based on a deterministic schedule, such as periodically vaccinating 
susceptible children based on school calendars. In reality, however, such policies are administered as a 
■ random process, while still possessing a mean period. Here, we consider the effect of randomly distributed 

intervention as disease control on large finite populations. We show explicitly how intervention control, 
based on mean period and treatment fraction, modulates the average extinction times as a function of 
population size and rate of infection spread. In particular, our results show an exponential improvement in 
extinction times even though the controls arc implemented using a random Poisson distribution. Finally, 
we discover those parameter regimes where random treatment yields an exponential improvement in 
extinction times over the application of strictly periodic intervention. The implication of our results is 
discussed in light of the availability of limited resources for control. 



i]n '. Introduction 

| Understanding the processes underlying disease extinction is an important problem in epidemic prediction 

. and control. Currently, total eradication of infectious disease is quite rare, but continues to be a major 

theme in public health. Temporary eradication, sometimes called fade out, tends to happen in local 
spatial regions, and may be followed by the reintroduction of the disease from other regions [U El E]. In 
the case of diseases that possess co-circulating strains such as influenza [4], or dengue fever which has 
up to four strains [5], extinction may occur in one or more strains while the others persist. Fade out 
or extinction may also occur within host, as is the case in Hepatitis C and HIV [BJ. Infectious disease 
transmission is also conjectured to be responsible for certain species extinction [7] [8]. Recently, large 
scale amphibian species have had major declines in population, which have been linked with the spread 
of disease [9]. 

One main reason that diseases go extinct is due to the stochasticity that is inherent to populations 
of finite size |10j . As a disease evolves in large finite population, there is the possibility of insufficient 
transmission for it to stay endemic. Therefore, in finite time, the number of infectious individuals can 
go to zero and the disease dies out [TTJ [12l [13] . Other mechanisms that enhance extinction include small 
populations and resource competition j!4j , as well as heterogeneity in host-vector models |15j . 

To properly model the random interactions occurring in populations, the study of disease extinction 
requires a stochastic modeling approach. As a result of the random interactions, time series analysis and 
epidemic models exhibit stochastic fluctuations [TBI H3 E] . The fluctuations may act as an effective 
force that drives the disease to vanish [20] . This force is composed of non-Markovian fluctuations which 
may overcome the instability of the extinct state. The non-Markovian nature is due to the pattern of the 
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noise necessary to drive the system out of equilibrium from the attracting state to the extinct state. It 
allows the escape from a quasi-stationary state by overcoming the force that keeps the system near the 
attractor. Therefore, if the fluctuations are weak on average, extinction is considered a rare event. 

We remark here that although escape has been considered for Langevin types of systems, the theory 
we present here is for discrete finite populations modeled as a master equation. In continuous systems, 
a rigorous theory of escape rates for systems driven by white Gaussian noise was developed by Freidlin 
and Wcntzell [21]. It was also found that the escape rates should display a number of universal features, 
including scaling behavior (221 Ei] , which has been confirmed by many experiments [24] EU [26] [27] . More 
recently, it was shown that the state of the system is coupled to a deterministic model of the noise shape 
[28) . In this setting, the optimal path is an unstable object, but may be associated with the dynamical 
systems idea of having maximum sensitive dependence to initial conditions [^5] . 

Vaccination and treatment programs are common methods used to speed up the extinction of a disease 
in a population |30j . In this paper, we aim to quantify how treatment or vaccination increase average 
extinction rates. We focus on a class of diseases with no immune response. Models with no immunity are 
suitable for many bacterial infections, such as meningitis, plague, and venereal diseases, as well as certain 
protozoan illnesses, such as malaria and sleeping sickness [31] . The mean field SIS model is typically 
analyzed for this scenario. If the disease is endemic, the disease state is attracting and the extinct, or 
disease free state, is unstable. 

Associated with the parameters for a particular disease is the basic reproduction number Rq, which 
defines on average how many new cases appear over one infectious period per infective pQ. When Rq < 1, 
the extinct state is attracting, but when Rq > 1, the extinct state is unstable. The bifurcation point at 
which the endemic state appears is at Rq=1. For example, it was approximated from the epidemiological 
data from England and Wales that the serogroup C meningococcal disease had Rq = 1-36, [32]. In Africa, 
some malaria Rq estimates are close to one, but others can be as high as 3,000 [33] . This variation is 
attributed to environmental temperature variations and mosquito biology [33]. Therefore, several groups 
have identified the applicability for methods to analyze extinction in finite populations near bifurcation 
points. (See the review in [35].) In both basic SIS [36] [37] and SIR [38] models, the mean times to 
extinction were analyzed as a function of Rq very close to one. The range of parameters here is assumed 
to model extremely slow disease propagation in large population limits. 

In general, little work has been done in analyzing stochastic models with random vaccination or 
treatment. Most vaccination schedules are designed as periodic, especially for childhood and seasonal 
diseases [39]. Each intervention typically has a prescribed (deterministic) schedule, or distribution, but 
the extinction event is still random. In [40] . a yearly vaccination pulse that was used as a control 
parameter was shown to prevent large amplitude disease outbreaks in a modified stochastic seasonal 
SEIR model. More recent work has considered the effect of deterministically imposed transitions on 
the rates of extinction [31]. Specifically, [3T] showed that a limited supply of vaccine can be optimally 
distributed to the susceptible population to change the rates of extinction exponentially in a range of low 

Rq. 

Thus, one of the main problems in understanding treatment and/or vaccination scheduling is that 
deterministic schedule models are not an accurate representation of the process. A more realistic scenario 
is that, on average, treatment scheduling has a mean period or cycle, but is itself a random process. In 
this paper, we study a randomly distributed treatment program of infected individuals. We are interested 
in evaluating treatment distributions by minimizing the mean time to extinction for the disease. Running 
simulations are computationally expensive and sensitive to population size. The theory presented here 
provides an alternate method to approximate the mean time to extinction. In our models, we identify 
conditions for which the escape rate theory applies and control strategies are effective. In particular, we 
derive explicit scaling functions of the exponent of the mean time to extinction in terms of Rq and mean 
treatment levels. We also identify the most effective treatment schedules. Then, we compare the theory 
against numerical simulations for verification. 
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Methods 

The stochastic SIS compartmental model tracks the number of individuals in a population of size N in one 
of two states: susceptible (Xi) or infected (X 2 ). In this model, we assume that the individuals become 
susceptible to the disease again upon recovery. The number of individuals in each state changes as birth, 
death, infection and recovery events occur. They are quantified by the following transition rates. New 
susceptible individuals are born at a rate [i, and both susceptible and infected individuals die at the rate 
5. In this model, we assume that the individuals recover from the disease without significant mortality. 
To keep the population constant over time, on average, we assume that the birth rate is equal to the 
death rate, so 6 = /i. If a susceptible comes in contact with an infected individual, the healthy individual 
may become infected. We use a mass action term with the contact rate j3 to describe the flow of newly 
infected individuals from the susceptible group. We assume infected individuals recover at rate k and 
immediately re-enter the susceptible group. In addition, we assume a treatment schedule that occurs at 
randomly chosen times with a frequency v times per year. Each time the treatment is applied, a fraction 
g of all infected individuals recover and flow back into the susceptible class. This assumes the treatment 
has 100% efficacy. To study the effect of treatments that are not as effective, a prefactor for g could be 
added to capture the smaller efficacy. That case is similar to studying a smaller value for g, which is 
included in the parameter range < g < 1 and therefore we do not study this issue separately. 

We now form the master equation that describes the time evolution of the stochastic system. The 
general theory of applying the WKB method to finite populations begins by assuming that the population 
of ./V individuals may be divided into n compartments. The number of people within each compartment 
is described by the n-dimcnsional state vector X = (Xi,X2, ■ ■ ■ ,X n ) with integer components. Let the 
random state transitions governing the dynamics be described by the transition rates W(X, r), with r 
representing the increment in the change of X, and r/N <C 1. Also, let the probability of finding the 
system in state X at time t be p(X, t). The general form of the master equation is 

MM = £[W(X - r; r)p(X - r, t) W(X; r)p(X, t)\. (1) 

r 

We assume that the system possess a single, strictly stationary solution, ^ = 0, that corresponds to the 
extinct state, where one or more of the n components of the vector X are equal to zero. 
In the large population limit (zero fluctuations), the mean field equations 

X = ^rV^(X;r) (2) 

r 

describe the time evolution of the system's mean values in a deterministic manner. 

The next step is to find a probability distribution, p(X, t), that solves the master equation. When the 
probability current at the extinct state is sufficiently small, there will exist a quasi-stationary probability 
distribution with a non-zero number of infected individuals that decays into the stationary solution over 
exponentially long times. The rate at which the extinction of infected individuals occurs may be calculated 
from the tail of the quasi-stationary distribution. It has been shown that a WKB approximation to the 
quasi-stationary distribution allows one to approximate the mean-time to extinction with high accuracy 
for N sufficiently large HH . This approximation consists in assuming that the desired distribution 
has, to leading order, the exponential form 

p(x,i) = Aexp(-7VS(x,<)), (3) 

where x = X/iV is the normalized state (e.g., in an epidemic model, the fraction of the population in the 
various compartments) and A is the normalization constant. The function S is given by an expansion of 
the form S = S + SJN + . . .. 
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Defining new scaled transition rates as w(x;r) = W(Nx;r)/N, inserting the ansatz of Eq. ^ into 
the master equation, Eq. ([1}, and keeping the first term only in the expansion for S yields a partial 
differential equation of the Hamilton- Jacobi form: 



at v ' 9x . 

In analogy to Hamiltonian mechanics, the functions H and S are called the Hamiltonian and the action, 
respectively. The derivative ^ is called the momentum conjugate to x and is denoted by p. With these 
approximations, the Hamiltonian function is 

H(x,p) =^w(x;r)[exp(p -r) - 1], (5) 

r 

and the characteristic equations for Eq. Q are precisely Hamilton's equations: 

x = d p H(x,p), p = -<9 x ff(x,p). (6) 
For a more detailed description of the WKB method and other applications, see [151 |4"4"1 155] . 

Model 1: Constrained SIS model with treatment 

In the first model, we make the following approximation in the SIS model to reduce the dimension of 
the problem. Assume the average population size is N and constrain the population size such that 
X\ + X 2 = N. Therefore, we can consider the dynamics of the constrained SIS model in terms of infected 
individuals, X 2 . Therefore, we need only to consider the following transition rates, which describe how 
individuals enter and leave the infected state: 

W(X 2 ;l) = fiX 2 (N-X 2 )/N, new infections; 
W(X 2 ; — l) = (p + k)X 2 , recovery and natural death; 

W(X 2 ] —gX 2 ) = v, treatment. 



Therefore, the master equation from Eq. (Q} for the constrained SIS stochastic process is 

(7) 



dp{X 2 .t) 
dt 



(ft + 7) ((X 2 + l)p(X 2 + 1, t) - X 2 p(X 2l i)) + u( K p{X 2 + gX 2 , t) - p(X 2 ,tfj 
+ £ ((X 2 - 1)(N - (X 2 - l))p{X 2 - 1, t) - X 2 (N - X 2 )p(X 2 ,t)) . 



Since the population variable in the master equation is integer- valued, the first argument of p(X 2 +gX 2 , t) 
must be rounded to an integer in a consistent way. Since treatment is assumed to reduce the infectious 
number of individuals, in what follows, we choose to keep the integer part [gX 2 ] of gX 2 (rounding down). 
Simplifying notation, we use gX 2 in the first argument of p(X 2 ,t) as shorthand for [3X2]. Note that 
for any particular realization of the master equation, the treatment ceases to have an effect whenever 
gX 2 < 1 (i.e., X 2 < 1/g) since for all those numbers of infecteds [gX 2 ] = 0. 
Following Eq. @, the associated deterministic model is 

X 2 = PX 2 {N - X 2 )/N -{p + K)X 2 -vgX 2 . (8) 

This system has two steady states: the extinct state, X 2 — 0, and the endemic state, X 2 = N(l — — ^S-), 
with Rq = j3 / {p, + k). As expected, more treatment (increasing g or v) decreases the number of infected 
individuals in the time-average. 

Equation always possesses as a solution a stationary distribution where the probability of observing 
zero infected individuals (X 2 = 0) is p(Q, t) = 1, which we identify as the extinct state. If Rq > 1 and N is 
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large enough, Eq. ([7]) will also possesses a quasi-stationary solution with an infected fraction fluctuating 
around an endemic state. Hence, if Rq > 1 the disease can spread through a population and is considered 
endemic. 

To simplify the variables for our model, we normalize the population so X2 = X2/N. Therefore, the 
Hamiltonian function for the SIS model is 

H(x 2 ,P2) = (3x 2 (l - x 2 ){e^ - 1) + Gu + n)x 2 {e-^ - 1) + L^-aN^ _ x) (9) 
The associated Hamiltonian system from Eq. © is 

X2 = I3x 2 (l-x 3 ) e P2 - (m + K)x 2 e-P* - vgx 2 e-3 Nx ^ , . . 

p 2 = -(3(1 -2x 2 ) (e P2 -l)-(p + K) (e~ P2 - 1) + vgp 2 e- gNx2P2 . [ ' 

For the Hamiltionian system, the endemic steady state is (x2 e ,P2 e ) = (1 — > 0) j while the stochastic 

die out state is (x2q, P2q) = (0,p*), with p* implicitly defined by 

vgp* = (3(e p * - 1) + (// + K )( e - p * - 1). (11) 

Note that the endemic state exists only if a; 2e > 0. In addition, the endemic state has zero momentum, 
which is consistent with our expectation that the probability distribution have a maximum at a^e & n d 
hence dS g^^ = P2e — 0. Since the variables X2 and P2 of the WKB approximation are not restricted 
to integer values, here the rounding of gX2 poses no problem. However, this means that in the WKB 
framework, the treatment pulses have an effect at arbitrarily low values of X2, in contrast to the master 
equation framework where, because of the rounding, treatment stops being applied whenever X2 < 1/5. 



Model 2: Full SIS model with treatment 

The second model is the unconstrained SIS treatment model in two-dimensions. We calculate S and / 
separately and allow the population fluctuation about N. If the fluctuations are small compared to N, 
the system will behave like the one-dimensional approximation. 

For the two-dimensional model, let the state vector be X = (Xi,X2) and the transition vector be 
r = (ri, 7-2). The changes in the susceptible and infected populations for a single transition are represented 
by the transition rates: 



W(X;(l,0y 


1 = Nix, 


birth of new susceptibles; 


W(X;(-1,0) 


\ = fxXt, 


natural death for susceptibles; 


W X; (0,-1) 


1 = px 2 , 


natural death for infectious; 


^(x;(i,-i); 


1 = nX 2l 


natural recovery; 


W(X;(-1,1); 


1 = pX x X 2 /N, 


new infections; 


W(X; {gX 2 , -gX 2 )) 


= 


treatment. 



Here, as in Model 1, the non-integer quantity gX2 is rounded down to [9X2]- We normalize the population 
so x = {x\,X2), with x\ = Xi/N and X2 = X2/N. Using the definition of the master equation, Eq. (JXJ) , 
and the ansatz for the probability density in Eq. ([3]), the Hamiltonian from Eq. ([5]) in normalized variables 
is 

H(x,p) = fi(e^ - 1) + (3 Xl X2{e-P 1+P2 - 1) + ra 2 (e pi - p2 - 1) + fj,xi(e- pi - 1) 
+fix 2 (e~P 2 + %( e 9 x * N Pi-'> 3 » N P* - 1). 

The associated Hamiltonian system from Eq. ((6|) is 



(12) 



±1 = neP 1 - (5 X!X2e- pi+ P 2 + kx 2 eP^'P 2 - pxie'P 1 + vgq 2 e 9X2N( - Pl -P 2 '> 

x 2 = (3 xix 2 e-P 1+P2 - KX 2 e Pl - p2 - px 2 e~ P2 - vgx2e gX2N ^-P 2 *> 

Pi = ~(3x2 (e~P 1+ P 2 - 1) -nie-P 1 - 1) 

p 2 = -(3x 1 (e-P^+P 2 - 1) - K(e Pl -P 2 - 1) - p{e- p2 - 1) - vg(p x - p 2 ) e ° x * N 



(13) 



G 



Note once more that in the WKB framework, the treatment pulses have an effect for arbitrarily small 
X2- For this Hamiltonian system, the endemic state is located at the point 

(^ 2 ,^ 2 )=(^ + f,l-^-f AO) (14) 

and the stochastic die out state is (a;i, X2,Pi,Pi) = (1,0, 0,p*) , with p* defined implicitly as in Eq. (fTTj) . 



Results 

We now use these Hamiltonian models to approximate the mean time to extinction. Topologically, the 
solution that describes an extinction event in the Hamiltonian system will connect the endemic state 
(x a ) and stochastic die out state (x s ). The connecting manifold is, in fact, the most probable path to 
extinction when the stochastic system starts initially at the endemic state [361 137] . This set of points 
is called the optimal path. Points on the path will also satisfy the Hamiltonian on the energy surface 
H(x,p) = since it is a solution to the time- independent version of Hamilton- Jacobi equation, Eq. ([6]). 

From the definition of the momentum, p(t) = the action up to the zeroth order of N along the 
optimal path can be found by evaluating 

Sopt = / Popt(t) ■ Xopt(t) dt. (15) 

Using this quantity, we approximate the mean time to extinction by 

r ext = Be NS °^ . (16) 

where B is a prcfactor that depends non-cxponentially on the system parameters and on the population 
size. An accurate approximation of the mean-time to extinction depends on obtaining B [45] . 

It is usually not a trivial task to identify the set of points that describe the optimal path. In some 
cases, it can be found analytically. One example is Model 1 with g = 0, since Eq. © has an explicit 
solution for p2 when constrained to H{x2,P2) = 0. An alternative approach is approximating the solution 
asymptotically. There are also several numerical approaches. One common method is to treat the system 
as a two point boundary value problem and solving using a shooting method |46j . In this paper, we use 
a generalized Newton's method that involves iterating an initial guess of the solution in the entire time 
domain [57j . Our initial guess must satisfy the property that the solution will stay asymptotically near the 
steady states except for a small, continuous, transition region between the two. This iterative procedure 
requires discretizing the model differential equations in time, using a second order approximation for the 
derivatives, and then solving the entire resulting system of nonlinear algebraic equations simultaneously. 

Equation (|16p holds if and only if a quasi-stationary distribution exists. This is the case if the time to 
extinction is exponentially long, i.e., NSopt 3> 1. Assuming that an endemic state does exist (x2e > 0); 
NS op t 3> 1 will be satisfied for N sufficiently large or, for fixed N, for an i?o sufficiently large and 
vg sufficiently small. The last conditions on the parameters mean that the disease should be highly 
transmissible and that the treatment should not be too intense. See the Supplementary Material, for a 
more detailed treatment on the necessary conditions for the quasi-stationary solution to exist. 



Model 1 

Because the Hamiltonian system for the constrainc 
to the action path simplifies to 




model is in two dimensions, the first approximation 
P2{x 2 )dx2, (17) 



7 



with P2 explicitly as a function of X2, evaluating the integral along the optimal path. The Hamiltonian 
function of Eq. ^ does not allow for an algebraic solution for ^2(^2) from the equation H{x2,P2) = 
that describes the path connecting the endemic state to the extinct state when g ^ 0. Therefore, this 
integral in Eq. (|17p must be approximated. 

For this model, an asymptotic approach can be used to approximate the action along the optimal 
path to extinction. We assume g is small, which implies small treatment pulses. We expand P2 in g and 
substitute this expression into the equation H(x2,P2) = 0. The resulting expansion is 

P2( X 2) = - HMl - (l - /3(1 _J g _ ( , + K) ) + °(9 2 )- (18) 

The first term in the expansion S — S + Si /N is given by Eq. ([T7|) . In (48] , the second term in the 
expansion of the action in powers of N is given as more complicated integral along the path. We expand 
the two integrals giving S and S\ in powers of g and evaluate them in closed form using computer algebra 
software. If we compare this asymptotic approximation to the numerical approximation for the action 
along the optimal path and evaluate Eq. (|16p . we see excellent agreement as shown in the example in 
Figure[T] In this example, we set the birth rate [i = 0.2 year -1 and recovery rate K = 100 year -1 . For the 
remainder of the paper, we will use these parameters in examples for both models. These values provide 
results that can be clearly visualized and easily reproduced. 

Note that while the action does not depend on the size of the population, to first order, the mean time 
to extinction does. The population size must be large enough to for the system to be quasi-stationary. 
Our model assumes that disease extinction is a rare event, which occurs in the tail of the distribution 
described by exp(— NS op t). Conversely, the peak of the distribution occurs at the endemic state. As Rq 
decreases to one, the distance between the endemic state and the disease free equilibrium decreases and 
the probability of the system having zero individuals in the infected state becomes significant. Therefore, 
the exponent must be large and negative, or cquivalently the action must be sufficiently large compared 
to the population so that S op t 3> 1/N. 

To quantify where the system is quasi-stationary, we use the numerical approximation of the optimal 
path to evaluate Eq. (fT7|) and construct a threshold graph to estimate the necessary parameters for the 
extinction to lie in the tail of the distribution. In Fig. [2] we show a threshold graph for the treatment 
model with frequency v = 4 year -1 . For a disease with j3 = 105, a population of 8, 000 is sufficiently to 
the right of the threshold for the exponent to be large for all < g < 0.4, meaning that the solution is 
quasi-stationary. 

The final step in finding the mean time to extinction is approximating the prcfactor in Eq. (j 16[) . 
Following the approach in [48] , we obtain^ 



B = I / 2itRo " 9 (19) 

(/?-( M+K )-^xi?o-i)y A^+«)in(i + ^r 1 > 

Note the dependence on the treatment parameters g and v. 

To quantify the accuracy of the approximation to the mean time to extinction in Eq. (|16|) . with S up 
to 0(N~ 1 ), we compare it to the average extinction time found by a Monte Carlo simulation as described 
in Gillespie |49j . In Figure [3] the graph shows this comparison over a range of treatment percentages 
(g) and frequencies (i>). The simulation uses a population of 10,000 and we averaged the results of 
2,000 realizations. As expected, the mean time to extinction decreases as the treatment percentage and 
frequency increase. Note the excellent agreement for small g, for which the asymptotic approximation 
was derived. 



1 We use Eq. (49) of I48| with A\ = 1 , which is the value that corresponds to our case. 
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Model 2 

The full SIS model lias a Hamiltonian system in four dimensions and asymptotic approximations of the 
optimal path and action are not tractable. Therefore, we rely on numerical approximations. We show 
projections of the locations of several paths in phase space in Figure |4j The arrow shows the direction 
of the path, connecting the endemic state to the extinct state. Notice how the paths decrease in length 
as we increase the control parameter g. Note that the action from integrating along those paths also 
decreases with g. Using these paths, we can approximate the action and the mean time to extinction for 
a treatment schedule. We compare the theory to data found by Monte Carlo simulation for small g, as 
shown in Figure [6] There is an exponential decrease in the mean time to extinction as we increase the 
treatment, agreeing with the theory. 

We also comment on the differences in the constrained and unconstrained SIS models with treatment. 
In Figure we compare the numerical approximations for actions of these models. We see that they 
agree for small Rq, but the action for the constrained model increases much faster as both Rq and g 
increase. This follows the result in [20], where the action for the constrained SIS model was shown to 
have an exponential scaling law for Rq ^> 1. The constrained treatment model also follows the exponential 
scaling. In addition, the theory can be used to avoid expensive simulations of long extinction times in 
large populations. The benefit of the full model is that it captures disease dynamics in a population with 
significant size fluctuations. The theory captures the rate of change in the mean time to extinction so 
that effectiveness in treatment schedules can be quantified. 

Discussion 

In this paper, we quantified how treatment enhances the extinction of epidemics using a stochastic, 
discrete-population framework. Specifically, we based our study on a general formulation of an SIS model 
with treatment that is applied randomly in a Poisson fashion, accounting for the limited amount of 
resources. We used a WKB approximation to the master equation of the stochastic process to calculate 
the average time to extinction starting from the endemic state, as a function of the transmissibility of 
the disease and the strength and frequency of the treatment. We compared the extinction times obtained 
analytically and numerically from the WKB approximation with the values obtained from Monte Carlo 
simulations. 

In addition, we explored the significance of the quasi-stationarity assumption that is fundamental to 
the WKB approximation. The existence of a quasi-stationary distribution peaked at the endemic point 
produces a meta-stable state in which the population fluctuates in a neighborhood around the same 
endemic point. In contrast, the extinct state lies in the exponentially small tail of the distribution. When 
a quasi-stationary distribution exists, the extinction of a disease is a rare event, i.e. the mean time to 
extinction is exponentially long. As we show in the Supplementary Material, the time to extinction is 
indeed exponentially long when the disease-free point lies in the tail of the distribution. The occurrence 
of extinction as a rare event means that the fluctuations exhibited by the random population dynamics 
are much smaller than an effective activation barrier. If the fluctuations are not small compared to the 
barrier, then the extinction events are not necessarily in the tail of the distribution, and hence not a rare 
event. 

Deterministic models of treatment and/or vaccination are not accurate representations of the process 
in finite population realizations. A more realistic description is that, on average, treatment scheduling 
has a mean period or cycle, but is itself a random process. To quantify the difference between the 
deterministic and the stochastic descriptions, we compared the mean time to extinction for a strictly 
periodic and a Poisson-distributcd treatment schedule obtained by averaging the Monte Carlo simulation 
results of many extinction events starting from the endemic state. We assume that a fraction g of the 
infected population is treated at a frequency of v times per year and immediately return to a susceptible 
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state. In Figure simulations support evidence that the random schedule had a faster mean time to 
extinction over the range of frequencies. The reason for this may be that when the system is close to the 
extinct state, there is a benefit to having a number of treatment pulses in a short window of time; such 
a series of frequent vaccination pulses are possible in the Poisson vaccination scheduling but not in the 
deterministic one. 

The treatment program that we implement in our model has two degrees of freedom: the frequency 
v and the fraction of infected individuals that are treated, g. On average, there are year) treatment 
pulses each year and at each one, a number Ngx2 of infected individuals are treated, where X2 is the 
infected fraction at the moment each treatment pulse occurs. Supposing that there are a fixed number 
of treatment doses Nvgx2e-{1 year)=constant that may be applied each year (here X2e is the fraction of 
the population that is infected at the endemic point). A natural question that arises is the following: 
Given a fixed number of total treatment doses, how are v and g chosen so that the time to disease 
extinction is minimized. In both of our SIS models, the fixed number of treatment doses translates into 
vg = constant. Monte Carlo simulations of Model 1 show that the mean time to extinction decreases 
uniformly as g increases, given a fixed vg quantity (Fig. [8]). The drop is particularly sharp for g — > 0. 
This appears to be a consequence of the rounding down of gX2 whenever a treatment pulse occurs (see 
Methods section). The treatment ceases to have an effect when there are less than 1/g infecteds; for very 
small g, the threshold 1/g is significant when compared to the number of infecteds at the endemic state. 
Thus, the treatment helps to bring the number of infected down to l/g, but not all the way to extinction. 
This issue does not appear if one instead chooses to round gX2 to the next-highest integer (results not 
shown). With this alternative method of rounding, the time to extinction actually has a sharp decrease 
as g — > 0. Monte Carlo simulations of Model 2 corroborate this finding. Thus, given a fixed number 
of resources, our stochastic simulations demonstrate that in order to eliminate infectious diseases, it is 
better to increase the pool of individuals reached by the treatment, rather than increase its frequency. 

In conclusion, we have described a method to quantify the effectiveness of a random treatment pro- 
gram. We find that increasing the magnitude and frequency of randomly scheduled treatments provide 
an exponential decrease in average extinction times. We have presented evidence that supports how 
larger campaigns applied less frequently are the most effective in facilitating disease eradication. Several 
assumptions in the model clarify the accuracy of the analytic approximation to the mean time to extinc- 
tion, but its exponential rate of decrease as we increase the intervention is consistent with simulations 
throughout our analysis as populations get very large. The techniques considered here can be easily gen- 
eralized to other diseases, such as those that include seasonality or population structure. Future work in 
this area could provide a more targeted control strategy that would be robust in fluctuating environments 
as well as more efficient and economical disease eradication. 
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Supplementary Material 

Quasi-stationarity 

In this section, we illustrate the central importance of the quasi-stationarity assumption for the accuracy 
of the WKB expression for the mean-extinction time. For this, we numerically compute solutions to the 
master equation for the constrained SIS model without treatment (Model 1 with g = 0), and determine 
when the WKB condition NS op t 3> 1 holds and when it does not by comparing the outcome with the 
WKB results. 

In the normalized constrained variables where x\-\-x% = 1 (i.e., no fluctuations in the total population), 
the Hamiltonian takes the form: 

H{x 2 ,p 2 ) = Px 2 {\ - x 2 )(e^ - 1) + {ii + n)x 2 {e~^ - 1). (20) 

The Hamiltonian system under the constraint H{x2-,P2) = has the steady states 

(x 2e , P 2e)= (i-^- (21) 

and (x 2 o,P2o) = (0, — lniio), which denote the endemic and extinct states, respectively. The system 
also possesses a third steady state (x2 m ,P2m) = (0,0). The attracting states {x2e,P2e) and (x2m,P2m) 
correspond to the zero fluctuation states that exist in the mean field equation (deterministic system). 
Since there is a nonzero probability current at the extinct state, (x2o,P2o) is a new state created by the 
noise in the system. We also note that for finite populations noise is not generally known due to the 
random interactions of individuals. However, in this model, the noise-free extinct state is accessible if 
the noise is known to be Gaussian. 

The most probable path from the endemic to the extinct point is the heteroclinic trajectory connecting 
the fixed point (x2 e ,P2e) with the fluctuational extinction point (x2o,p 2 o)- The optimal path is given by 
^2(2^2) = — hi(i?o(l — x 2.))- Therefore, the action from the endemic state to a point X2 along the optimal 
path up to the zeroth order of N is 

S{x 2 )= -ln(Ro(l-X2))dx 2 . (22) 

In particular, the action from the endemic state to the extinct state along the optimal path is 

S opt = S(0) = -l + la(Ro) + ^-- (23) 

Since Rq > 1, for sufficiently large N we have NS op t 3> 1, ensuring the extinct point lies in the tail of 
the probability distribution, where its value p(0) = Ae~ NS ° pt is exponentially small. 

In Figure O the WKB approximations to the quasi-stationary distribution are shown for the cases of 
R = 2 and Rq = 1.1. The value of the distribution at zero shows that the extinction point is definitely 
not in the tail of the distribution for Rq = 1.1 and hence does not constitute a rare event. In contrast, 
for i?o = 2, the extinct state is in the tail of the distribution and hence we expect the WKB result to be 
accurate. 

Because the disease free equilibrium is always an absorbing boundary in the one-dimensional case, 
the system decays to the extinct state in the long term. If a quasi-stationary distribution exists, the 
complete decay happens for exponentially long times; otherwise, it occurs on a much shorter time-scale. 
To illustrate this phenomenon, we show the numerical solution of the master equation in Eq. (JT)) over 
time in Figure [TO] The initial probability distribution at t = is set to the WKB approximation of the 
SIS probability distributions using Eq. (f2"2"j) . The absorption into the disease- free state is apparent in the 
R = 1.1 case, but completely imperceptible for Rq = 2 over the time-scale shown. 
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Figure Legends 




Figure 1. Comparing quantitative approximations of the action. For Model 1, plot of the 
numerical approximation of the action (dashed curve) and the asymptotic approximation (solid curve) 
as a function of the treatment, g. In this example, we use the parameters j3 — 105 year -1 and N = 8000 
people. As expected, the best agreement is for small g. 




Figure 2. Checking the threshold for quasi-stationarity. A plot of the log of the action for 
Model 1 as we vary treatment, g, vs. the contact rate, j3. In this case, the treatment frequency is v = 4 
year -1 . The red line indicates the threshold 1/N for N = 8, 000 people. The action for (3 = 105 year -1 
and N = 8, 000 people is greater than the threshold as we vary g. Therefore, we conclude these 
parameters would be sufficient for the model to exhibit extinction in the tail of the probability. As we 
increase the treatment fraction g for /3 = 105 year -1 , the action decreases towards the threshold line 
and a larger population would be necessary. 
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Figure 3. The effectiveness of various treatment combinations for Model 1. A plot of the 
fraction of infected treated, g, vs. the mean time to disease extinction, r ext years, for different 
treatment frequencies, v year -1 . The results for the Monte Carlo simulations are averaged over 2,000 
realizations and plotted as symbols. The curves of the same color show the approximation of the mean 
time to extinction by finding the action. The parameters are N = 8, 000 people and (3 = 105 year -1 . 
Note the exponential decrease in the mean time to extinction as the treatment fraction is increased. 




Figure 4. Graphs of the numerical approximations of optimal paths for Model 2. In (a), we 
show the (xi,X2,Pi) projection. In (b), we show the (xi, £2,.P2) projection. The parameters are v = 4 
year -1 , (3 = 105 and g = 0,0.15,0.3 for these examples. The arrow shows the direction of the path, 
connecting the endemic state to the extinct state. Below the curves are the projections of the x\ and X2 
onto the respective momenta planes. 



16 




Figure 5. A comparison of the action approximations for Model 1 and Model 2. This plot 
shows the quantitative difference in the action approximation for Model 1 (blue) and Model 2 (red) as 
we vary Rq and g. In this example, we use parameters /3 = 105 year -1 , v = 4 year -1 , and N = 120, 000. 




Figure 6. The effectiveness of various treatment combinations for Model 2. A plot of the 
fraction of infected vaccinated during each treatment, g, vs. the mean time to extinction, T ext , for 
different treatment frequencies, v year -1 . The average of 2,000 Monte Carlo simulations are shown with 
symbols. The curves of the same color show the numerical approximation of T ext using the action S and 
a constant prefactor. For the parameters, we use N = 12, 000 people and /3 = 105 year -1 . 
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Figure 7. A comparison of periodic and random treatment effectiveness. For Model 1, a plot 
of the fraction of infected vaccinated during each treatment vs. the Monte Carlo simulated mean time 
to disease extinction for random (points connected by dotted lines) and periodic (symbols) treatment 
schedules. Results arc shown for treatment frequencies, v =2, 4, 8, and 12 year -1 averaged over 2,000 
realizations. The parameters are N = 8,000 people and f3 = 105 year -1 . Note that the random 
treatment schedule has average extinction times consistently lower than the periodic treatment schedule. 
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Figure 8. The effectiveness of various treatment combinations for a fixed treatment 
supply. For Model 1, a plot of the mean time to extinction as we vary g for a fixed treatment supply 
gv =constant. The symbols represent the Monte Carlo simulation results for gv =0.1, 0.2, and 0.3 
averaged over 10,000 realizations. The curves represent the direct numerical solution of the associated 
master equation. The parameters are N = 8, 000 people and {3 = 105 year -1 . Note the decrease in 
mean time to extinction as the g increases. 
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Figure 9. Quasi-stationarity depicted through probability distributions. Graphs of the WKB 
approximation of the SIS probability distributions using Eq. (f22j) for N = 200. We show the case of 
Rq = 2, for which extinction is in the tail of the distribution. Conversely, extinction has a significant 
probability in the case of Rq = 1.1. Note the height of the curve for X2 = 0. The dotted vertical lines 
show the location of the endemic state in Xi for each case. 




Figure 10. The drift of probability distributions for systems without quasi-stationarity. A 

plot of the solution of the ID master equation in Eq. with g = over time using the distribution 
from the WKB approximation, Eq. (|22|) . as initial condition. For Rq = 2 (left), the extinct state lies in 
the tail of the distribution and a quasi-stationary distribution exists. Extinction occurs only over 
exponentially long times. For Rq = 1.1 (right) the endemic state is close to the absorbing boundary and 
extinction is not a rare event. The absorption of this distribution into the boundary is apparent. 



